#include "evaluate.h"
#include <cstring>
#include <vector>
#include <algorithm>
#include <numeric>
#include <utility>
#include <cmath>

using namespace std;

void GetEvaluation(
    double* x, int* y, int n, int nCluster, double* centers, double* size,
    int* nScore, double** pScore, char*** pMethod) {
    static const int MAXM = 12;
    static const int MAXS = 80;
//    static const int MAXN = 6e5;
    static const int MAXK = 10;
    static const char METHOD[MAXM][MAXS] = {
        "(-内聚) Inertia / SSE 簇内误方差",
        "(-内聚) Compactness 紧密性",
        "(+耦合) Separation 间隔性",
        "(-内聚耦合) Davies-Bouldin Index 戴维森堡丁指数",
        "(+内聚耦合) Dunn Validity Index 邓恩指数",
        "(+内聚耦合) Silhouette coefficient 轮廓系数",
        "(+内聚耦合) Calinski-Harabaz / Variance Ratio Criterion 方差比标准",
        "(-均匀) Variance 方差",
        "(-均匀) Sample standard deviation 样本标准偏差",
        "(-均匀) Coefficient of variance 变异系数",
        "(+均匀) Information Entropy 信息熵",
        "(+均匀) Gini Index 基尼指数"
    };

    *nScore = MAXM;
    *pMethod = new char* [MAXM];
    char** method = *pMethod;
    for (int i = 0; i < MAXM; i++) {
        method[i] = new char[MAXS];
        strcpy(method[i], METHOD[i]);
    }
    *pScore = new double[MAXM];
    double* score = *pScore;

    static vector<double> cluster[MAXK];
    double sumX=0;
    for (int i = 0; i < nCluster; i++) cluster[i].clear();
    for(int i=0; i<n; i++) {
        cluster[y[i]].push_back(x[i]);
        sumX+=x[i];
    }
    static vector<double> sum[MAXK];
    static vector<double> sum2[MAXK];
    static double innerS[MAXK];
    static double innerS2[MAXK];
    static double innerUS[MAXK];
//    static double innerUS2[MAXK];
    static double sumCenter[MAXK+1];
    sumCenter[0]=0;
    for(int i=0; i<nCluster; i++) {
        vector<double>& c = cluster[i];
        sort(c.begin(),c.end());
        sum[i].resize(c.size()+1);
        sum2[i].resize(c.size()+1);
        sum[i][0]=0;
        sum2[i][0]=0;
        for (size_t j = 0; j < c.size(); j++) {
            sum[i][j+1]=sum[i][j]+c[j];
            sum2[i][j+1]=sum2[i][j]+c[j]*c[j];
        }
        int ub=upper_bound(c.begin(),c.end(),centers[i])-c.begin();
        innerS[i]=sum[i][size[i]]-sum[i][ub]-centers[i]*(size[i]-ub)
                  +centers[i]*ub-(sum[i][ub]-sum[i][0]);
        innerUS[i]=innerS[i]/size[i];
        double a=centers[i],a2=a*a*size[i];
        double b=sum[i][size[i]]-sum[i][0],b2=sum2[i][size[i]]-sum2[i][0];
        innerS2[i]=a2-2*a*b+b2;
//        innerUS2[i]=innerS2[i]/size[i];
        sumCenter[i+1]=sumCenter[i]+centers[i];
    }

    // 0: (-内聚) Inertia / SSE 簇内误方差
    score[0]=0;
    for(int i=0; i<nCluster; i++) score[0]+=innerS2[i];
    // 1: (-内聚) Compactness 紧密性
    score[1]=0;
    for(int i=0; i<nCluster; i++) score[1]+=innerS[i]/size[i];
    score[1]/=nCluster;
    // 2: (+耦合) Separation 间隔性
    score[2]=0;
    for(int i=1; i<=nCluster; i++)
        score[2]+=sumCenter[nCluster]-sumCenter[i]-centers[i-1]*(nCluster-i);
    score[2]/=nCluster*(nCluster-1)/2;
    // 3: (-内聚耦合) Davies-Bouldin Index 戴维森堡丁指数
    score[3]=0;
    for(int i=0; i<nCluster; i++) {
        double up=0;
        for(int j=0; j<nCluster; j++) if(i!=j)
                up=max(up,(innerUS[i]+innerUS[j])/abs(centers[i]-centers[j]));
        score[3]+=up;
    }
    score[3]/=nCluster;
    // 4: (+内聚耦合) Dunn Validity Index 邓恩指数
    double minIntraPair=cluster[1][0]-cluster[0][size[0]-1];
    double maxInnerPair=cluster[0][size[0]-1]-cluster[0][0];
    for(int i=1; i<nCluster; i++) {
        vector<double> &c=cluster[i];
        minIntraPair=min(minIntraPair,c[0]-cluster[i-1][size[i-1]-1]);
        maxInnerPair=max(maxInnerPair,c[size[i]-1]-c[0]);
    }
    score[4]=minIntraPair/maxInnerPair;
    // 5: (+内聚耦合) Silhouette coefficient 轮廓系数
    score[5]=0;
    double xi,ai,bi,si;
    for(int c=0; c<nCluster; c++) {
        for(int i=0; i<size[c]; i++) {
            xi=cluster[c][i];
            ai=(sum[c][size[c]]-sum[c][i+1]-xi*(size[c]-1-i)
                +xi*i-(sum[c][i]-sum[c][0]))/(size[c]-1);
            if(c==0) bi=centers[c+1]-xi;
            else if(c==nCluster-1) bi=xi-centers[c-1];
            else bi=min(centers[c+1]-xi,xi-centers[c-1]);
            si=(bi-ai)/max(ai,bi);
            score[5]+=si;
        }
    }
    score[5]/=n;
    // 6: (+内聚耦合) Calinski-Harabaz / Variance Ratio Criterion 方差比标准
    double w=accumulate(innerS2,innerS2+nCluster,0.0);
    double ux=sumX/n;
    double b=0;
    for(int i=0; i<nCluster; i++)
        b+=size[i]*pow(abs(centers[i]-ux),2);
    score[6]=b/w*(n-nCluster)/(nCluster-1);
    // 7: (-均匀) Variance 方差
    double sumSize=accumulate(size,size+nCluster,0.0);
    double uSize=sumSize/nCluster;
    score[7]=0;
    for(int i=0; i<nCluster; i++) score[7]+=pow(size[i]-uSize,2);
    score[7]/=nCluster;
    // 8: (-均匀) Sample standard deviation 样本标准偏差
    score[8]=sqrt(score[7]*nCluster/(nCluster-1));
    // 9: (-均匀) Coefficient of variance 变异系数
    score[9]=sqrt(score[7])/uSize;
    // 10: (+均匀) Information Entropy 信息熵
    score[10]=0;
    for(int i=0; i<nCluster; i++)
        score[10]-=size[i]/n*(log2(size[i])-log2(n));
    // 11: (+均匀) Gini Index 基尼指数
    score[11]=1;
    for(int i=0; i<nCluster; i++) score[11]-=pow(size[i]/n,2);
}
